Antje Hofmann: Learning Syntactic Relations with the Mole Task (2)

RePsychLing in SMLP2022

Author

Reinhold Kliegl

1 Background

This version removes the outlier subject “p10” found in the version 1.

Children’s (age: 4-8 years)reaction times in a task teaching them syntactic relations.

1.1 Overview

  • Original analysis is by Antje Hofmann.
  • MixedModels.jl version
  • Addition of new chunks illustrate
    • selection of parsimonious LMM using random-effects PCA
    • plotting conditional means
    • illustration of borrowing strength

2 Readme

2.1 Variables

  • Subj: Participant ID (renamed from ID; random factor)
  • Item: Word ID (random factor)
  • age: 4 - 8 years
  • Block (within-Subj/within-Item):
    • 1st Learning
    • 2nd Learning
    • Disruption
    • Recovery
  • Target(renamend fom targetness)
    • non-syllable target
    • syllable target
  • rt: response time

3 Setup

3.1 Packages

First attach the MixedModels.jl package and other packages for plotting. The CairoMakie.jl package allows the Makie graphics system [@Danisch2021] to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.

Code
using AlgebraOfGraphics
using Arrow
using CairoMakie       # graphics back-end
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros  # simplified dplyr-like data wrangling
using MixedModels
using MixedModelsMakie # diagnostic plots
using ProgressMeter
using Random           # random number generators
using RCall            # call R from Julia
using StatsModels

using AlgebraOfGraphics: boxplot
using AlgebraOfGraphics: density

using MixedModelsMakie: qqnorm
using MixedModelsMakie: ridgeplot
using MixedModelsMakie: scatter
using MixedModelsMakie: caterpillar

ProgressMeter.ijulia_behavior(:clear);
CairoMakie.activate!(; type="svg");
  • The data are available as an arrow file.
  • Most of preprocessing was done with R in RStudio (see Hofmann_Maulwurf.Rmd).
  • Order of factor levels should be checked.
Code
dat = DataFrame(Arrow.Table("./data/Hofmann_Maulwurf_rt.arrow"))
transform!(dat,
    :Target => categorical => :Target,
    :Block  => categorical => :Block,
    :age => (x -> x .- 6) => :a1, # center age at six years
    :rt => (x -> log.(x)) => :lrt)
describe(dat)
dat = filter(row -> row.Subj != "p10", dat)

6,257 rows × 8 columns (omitted printing of 1 columns)

SubjItemageBlockTargetrta1
String?String?Float64?Cat…?Cat…?Float64?Float64
1p0102SG.ogg7.731st LearningNon-target syllable3285.01.73
2p0121SG.ogg7.731st LearningNon-target syllable2873.01.73
3p0120SG.ogg7.731st LearningNon-target syllable2169.01.73
4p0113SG.ogg7.731st LearningNon-target syllable3133.01.73
5p0132SG.ogg7.731st LearningNon-target syllable5793.01.73
6p0105SG.ogg7.731st LearningNon-target syllable2064.01.73
7p0127SG.ogg7.731st LearningNon-target syllable4546.01.73
8p0118SG.ogg7.731st LearningNon-target syllable2949.01.73
9p0123SG.ogg7.731st LearningNon-target syllable2957.01.73
10p0107SG.ogg7.731st LearningNon-target syllable2508.01.73
11p0129SG.ogg7.731st LearningNon-target syllable2232.01.73
12p0119SG.ogg7.731st LearningNon-target syllable2517.01.73
13p0114SG.ogg7.731st LearningNon-target syllable2234.01.73
14p0116SG.ogg7.731st LearningNon-target syllable4272.01.73
15p0115SG.ogg7.731st LearningNon-target syllable3389.01.73
16p0109SG.ogg7.731st LearningNon-target syllable2510.01.73
17p0125SG.ogg7.731st LearningNon-target syllable5079.01.73
18p0112SG.ogg7.731st LearningNon-target syllable2622.01.73
19p0123PA.ogg7.731st LearningTarget syllable2161.01.73
20p0132PA.ogg7.731st LearningTarget syllable2141.01.73
21p0119PA.ogg7.731st LearningTarget syllable2311.01.73
22p0113PA.ogg7.731st LearningTarget syllable3049.01.73
23p0104PA.ogg7.731st LearningTarget syllable2175.01.73
24p0125PA.ogg7.731st LearningTarget syllable1750.01.73
25p0120PA.ogg7.731st LearningTarget syllable2158.01.73
26p0114PA.ogg7.731st LearningTarget syllable3274.01.73
27p0111PA.ogg7.731st LearningTarget syllable2434.01.73
28p0130PA.ogg7.731st LearningTarget syllable1859.01.73
29p0115PA.ogg7.731st LearningTarget syllable2351.01.73
30p0116PA.ogg7.731st LearningTarget syllable1853.01.73
  • Centering age at six years yields an interpretable GM
  • Factor levels can also be set when contrasts are defined (see below).
  • BoxCox check showed that reaction time rt [ms] should be transformed to speed [1/s] = [Hz]
  • Indicator variables for Target and Block generated in R.

4 LMM analysis

4.1 Contrasts

Code
contrasts = merge(
      Dict(:Target => EffectsCoding()),
      Dict(:Block => SeqDiffCoding()),
      Dict(nm => Grouping() for nm in (:Subj, :Item))
   );

4.2 Varying only intercepts LMM m_voi

Code
f_voi1 = @formula(lrt ~  1 + Block * Target * a1 + (1 | Subj) + (1 | Item));
m_voi1 = fit(MixedModel, f_voi1, dat; contrasts)
Minimizing 48    Time: 0:00:01 (24.56 ms/it)
Est. SE z p σ_Item σ_Subj
(Intercept) 7.7048 0.0358 215.08 <1e-99 0.0255 0.2508
Block: 2nd Learning -0.0810 0.0102 -7.95 <1e-14
Block: Disruption 0.0518 0.0124 4.18 <1e-04
Block: Recovery -0.0509 0.0123 -4.14 <1e-04
Target: Target syllable -0.0121 0.0040 -3.05 0.0023
a1 -0.1054 0.0314 -3.36 0.0008
Block: 2nd Learning & Target: Target syllable -0.0003 0.0102 -0.03 0.9746
Block: Disruption & Target: Target syllable -0.0383 0.0122 -3.15 0.0016
Block: Recovery & Target: Target syllable 0.0102 0.0121 0.85 0.3959
Block: 2nd Learning & a1 0.0237 0.0089 2.65 0.0081
Block: Disruption & a1 0.0211 0.0106 1.99 0.0470
Block: Recovery & a1 -0.0127 0.0106 -1.20 0.2305
Target: Target syllable & a1 -0.0040 0.0035 -1.15 0.2494
Block: 2nd Learning & Target: Target syllable & a1 -0.0165 0.0089 -1.85 0.0641
Block: Disruption & Target: Target syllable & a1 0.0234 0.0106 2.21 0.0274
Block: Recovery & Target: Target syllable & a1 0.0005 0.0105 0.04 0.9656
Residual 0.2919

4.3 Extract indicator variables

Code
X = modelmatrix(m_voi1)
dat.trng = X[:,2];
dat.drpt = X[:,3];
dat.rcvr = X[:,4];
dat.trgt = X[:,5];
describe(dat)

12 rows × 7 columns

variablemeanminmedianmaxnmissingeltype
SymbolUnion…AnyUnion…AnyInt64Type
1Subjp01p850Union{Missing, String}
2Item01PA.ogg32SG.ogg0Union{Missing, String}
3age6.231143.76.318.040Union{Missing, Float64}
4Block1st LearningRecovery0Union{Missing, CategoricalValue{String, UInt32}}
5TargetNon-target syllableTarget syllable0Union{Missing, CategoricalValue{String, UInt32}}
6rt2343.97228.02105.08389.00Union{Missing, Float64}
7a10.231143-2.30.312.040Float64
8lrt7.679665.429357.652079.034680Float64
9trng-0.0140243-0.750.250.250Float64
10drpt-0.0518619-0.5-0.50.50Float64
11rcvr0.0530206-0.25-0.250.750Float64
12trgt0.00175803-1.01.01.00Float64

Switch to indicator variables and refit the model.

Code
f_voi2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1 + (1 | Subj) + (1 | Item));
m_voi2 = fit(MixedModel, f_voi2, dat; contrasts)
Est. SE z p σ_Item σ_Subj
(Intercept) 7.7048 0.0358 215.08 <1e-99 0.0255 0.2508
trng -0.0810 0.0102 -7.95 <1e-14
drpt 0.0518 0.0124 4.18 <1e-04
rcvr -0.0509 0.0123 -4.14 <1e-04
trgt -0.0121 0.0040 -3.05 0.0023
a1 -0.1054 0.0314 -3.36 0.0008
trng & trgt -0.0003 0.0102 -0.03 0.9746
drpt & trgt -0.0383 0.0122 -3.15 0.0016
rcvr & trgt 0.0102 0.0121 0.85 0.3959
trng & a1 0.0237 0.0089 2.65 0.0081
drpt & a1 0.0211 0.0106 1.99 0.0470
rcvr & a1 -0.0127 0.0106 -1.20 0.2305
trgt & a1 -0.0040 0.0035 -1.15 0.2494
trng & trgt & a1 -0.0165 0.0089 -1.85 0.0641
drpt & trgt & a1 0.0234 0.0106 2.21 0.0274
rcvr & trgt & a1 0.0005 0.0105 0.04 0.9656
Residual 0.2919

They are equivalent.

4.4 A zero-correlation parameter LMM m_zcp

Code
f_zcp1 = @formula(lrt ~  1 +  Block * Target * a1 + zerocorr(1 + Block * Target | Subj) + (1 + a1 | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat; contrasts);

show(issingular(m_zcp1))
VarCorr(m_zcp1)
Minimizing 375   Time: 0:00:00 ( 1.46 ms/it)
false
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0668442 0.2585425
Block: 2nd Learning 0.0305112 0.1746747 .
Block: Disruption 0.0165614 0.1286909 . .
Block: Recovery 0.0135351 0.1163403 . . .
Target: Target syllable 0.0002921 0.0170918 . . . .
Block: 2nd Learning & Target: Target syllable 0.0003617 0.0190193 . . . . .
Block: Disruption & Target: Target syllable 0.0023922 0.0489103 . . . . . .
Block: Recovery & Target: Target syllable 0.0006247 0.0249938 . . . . . . .
Item (Intercept) 0.0009965 0.0315682
a1 0.0005789 0.0240613 -0.77
Residual 0.0724795 0.2692202

Again, check the equivalence.

Code
f_zcp2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1  +
                zerocorr(1 + (trng+drpt+rcvr) * trgt | Subj) + (1 + a1 | Item));
m_zcp2 = fit(MixedModel, f_zcp2, dat; contrasts);

show(issingular(m_zcp2))
VarCorr(m_zcp2)
Minimizing 375   Time: 0:00:00 ( 0.78 ms/it)
false
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0668442 0.2585425
trng 0.0305112 0.1746747 .
drpt 0.0165614 0.1286909 . .
rcvr 0.0135351 0.1163403 . . .
trgt 0.0002921 0.0170918 . . . .
trng & trgt 0.0003617 0.0190193 . . . . .
drpt & trgt 0.0023922 0.0489103 . . . . . .
rcvr & trgt 0.0006247 0.0249938 . . . . . . .
Item (Intercept) 0.0009965 0.0315682
a1 0.0005789 0.0240613 -0.77
Residual 0.0724795 0.2692202

4.5 A complex parameter LMM m_cpx

Code
m_cpx = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt * (trng+drpt+rcvr) | Subj) + (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_cpx))    # not ok
show(m_cpx.PCA[:Subj])     # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_cpx))
VarCorr(m_cpx)
Minimizing 1640      Time: 0:00:01 ( 0.64 ms/it)
  objective:  1954.9673786896176
true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .     .
 trgt         -0.34   1.0     .      .      .      .      .     .
 trng          0.23   0.13   1.0     .      .      .      .     .
 drpt         -0.01  -0.08   0.3    1.0     .      .      .     .
 rcvr         -0.35   0.01  -0.13  -0.45   1.0     .      .     .
 trgt & trng   0.06  -0.52  -0.49   0.01   0.02   1.0     .     .
 trgt & drpt  -0.16   0.06   0.44   0.42  -0.38  -0.64   1.0    .
 trgt & rcvr   0.54   0.26  -0.36  -0.16  -0.09   0.13  -0.6   1.0

Normalized cumulative variances:
[0.33, 0.5515, 0.7424, 0.8685, 0.9466, 0.9934, 1.0, 1.0]

Component loadings
 
               PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8
 (Intercept)  -0.12   0.6    0.28  -0.44   0.08   0.11   0.34  -0.47
 trgt          0.14  -0.33   0.59   0.39  -0.21  -0.29   0.39  -0.29
 trng          0.43   0.13   0.14  -0.48  -0.51  -0.39  -0.07   0.36
 drpt          0.31   0.36  -0.2    0.43  -0.57   0.43  -0.06  -0.19
 rcvr         -0.21  -0.54  -0.08  -0.41  -0.38   0.51   0.3   -0.02
 trgt & trng  -0.43   0.2   -0.45   0.18  -0.24  -0.39   0.55   0.2
 trgt & drpt   0.57   0.03  -0.09   0.04   0.4    0.25   0.57   0.34
 trgt & rcvr  -0.37   0.24   0.55   0.18  -0.08   0.31  -0.02   0.6
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         56  1954.9674  52.1807      28  0.0037
──────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0676775 0.2601489
trgt 0.0003092 0.0175844 -0.34
trng 0.0288450 0.1698380 +0.23 +0.13
drpt 0.0196619 0.1402209 -0.01 -0.08 +0.30
rcvr 0.0180750 0.1344433 -0.35 +0.01 -0.13 -0.45
trgt & trng 0.0015803 0.0397531 +0.06 -0.52 -0.49 +0.01 +0.02
trgt & drpt 0.0065055 0.0806569 -0.16 +0.06 +0.44 +0.42 -0.38 -0.64
trgt & rcvr 0.0030635 0.0553486 +0.54 +0.26 -0.36 -0.16 -0.09 +0.13 -0.60
Item (Intercept) 0.0008437 0.0290467
a1 0.0005576 0.0236127 -0.87
Residual 0.0720702 0.2684589

The deviance improves, but we end up with an overparameterized LMM.

4.6 A parsimonious parameter LMM m_prm

We remove one of the VC for trgt * trng contrast interaction, that is one of the three interaction terms.

Code
m_prm1 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt *(drpt + rcvr) + trng  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm1))  # ok
show(m_prm1.PCA[:Subj])   # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm1, m_cpx))
VarCorr(m_prm1)
Minimizing 2205      Time: 0:00:01 ( 0.63 ms/it)
true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .
 trgt         -0.33   1.0     .      .      .      .      .
 drpt         -0.01  -0.11   1.0     .      .      .      .
 rcvr         -0.34   0.01  -0.45   1.0     .      .      .
 trgt & drpt  -0.18  -0.09   0.5   -0.43   1.0     .      .
 trgt & rcvr   0.54   0.28  -0.16  -0.09  -0.67   1.0     .
 trng          0.21   0.08   0.3   -0.12   0.4   -0.38   1.0

Normalized cumulative variances:
[0.3473, 0.5941, 0.7601, 0.8926, 0.9713, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.1    0.68  -0.24   0.32  -0.11   0.3   -0.52
 trgt         -0.12  -0.16   0.87   0.19  -0.1    0.13  -0.37
 drpt          0.44   0.23   0.17  -0.28   0.78   0.02  -0.19
 rcvr         -0.29  -0.5   -0.26   0.34   0.46   0.51  -0.07
 trgt & drpt   0.58  -0.05   0.05  -0.11  -0.33   0.7    0.22
 trgt & rcvr  -0.46   0.45   0.28  -0.07   0.19   0.29   0.62
 trng          0.38   0.1    0.08   0.8    0.12  -0.24   0.34
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + drpt + rcvr + trgt & drpt + trgt & rcvr + trng | Subj) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         48  1962.8756  44.2725      20  0.0014
[3]         56  1954.9674   7.9082       8  0.4425
──────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0678202 0.2604231
trgt 0.0003055 0.0174778 -0.33
drpt 0.0195723 0.1399012 -0.01 -0.11
rcvr 0.0180836 0.1344754 -0.34 +0.01 -0.45
trgt & drpt 0.0049548 0.0703907 -0.18 -0.09 +0.50 -0.43
trgt & rcvr 0.0029429 0.0542485 +0.54 +0.28 -0.16 -0.09 -0.67
trng 0.0297813 0.1725726 +0.21 +0.08 +0.30 -0.12 +0.40 -0.38
Item (Intercept) 0.0008334 0.0288684
a1 0.0005466 0.0233798 -0.87
Residual 0.0722829 0.2688548

We don’t lose goodness of fit, but are still overparameterized. We remove CPs for Target.

Code
m_prm2 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm2))  #  ok
show(m_prm2.PCA[:Subj])   #  ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm2,  m_cpx))
VarCorr(m_prm2)
Minimizing 1064      Time: 0:00:00 ( 0.58 ms/it)
false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt         -0.01   1.0     .      .      .      .     .
 rcvr         -0.34  -0.45   1.0     .      .      .     .
 trng          0.21   0.3   -0.12   1.0     .      .     .
 drpt & trgt  -0.12   0.53  -0.44   0.42   1.0     .     .
 rcvr & trgt   0.5   -0.2   -0.09  -0.44  -0.67   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3523, 0.5915, 0.7343, 0.869, 0.9457, 0.9836, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.08   0.65  0.0   -0.44   0.14  -0.52  -0.3
 drpt          0.45   0.22  0.0    0.34  -0.76  -0.17  -0.18
 rcvr         -0.3   -0.54  0.0   -0.36  -0.41  -0.54   0.17
 trng          0.41   0.07  0.0   -0.73  -0.22   0.44   0.23
 drpt & trgt   0.58  -0.03  0.0    0.11   0.35  -0.46   0.56
 rcvr & trgt  -0.46   0.48  0.0    0.14  -0.25   0.09   0.69
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         42  1969.2274  37.9206      14  0.0005
[3]         56  1954.9674  14.2600      14  0.4305
──────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0678705 0.2605197
drpt 0.0195960 0.1399859 -0.01
rcvr 0.0181233 0.1346227 -0.34 -0.45
trng 0.0298994 0.1729145 +0.21 +0.30 -0.12
drpt & trgt 0.0048302 0.0694998 -0.12 +0.53 -0.44 +0.42
rcvr & trgt 0.0024143 0.0491356 +0.50 -0.20 -0.09 -0.44 -0.67
trgt 0.0003211 0.0179201 . . . . . .
Item (Intercept) 0.0008544 0.0292300
a1 0.0005640 0.0237488 -0.83
Residual 0.0723240 0.2689312

Check the Item-related CP. It is very large, might be spurious.

Code
m_prm3 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                  zerocorr(1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm3))  # ok
show(m_prm3.PCA[:Subj])   # ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm3, m_prm2, m_cpx))
VarCorr(m_prm3)
Minimizing 1133      Time: 0:00:00 ( 0.55 ms/it)
false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt         -0.01   1.0     .      .      .      .     .
 rcvr         -0.34  -0.45   1.0     .      .      .     .
 trng          0.21   0.31  -0.11   1.0     .      .     .
 drpt & trgt  -0.13   0.53  -0.44   0.43   1.0     .     .
 rcvr & trgt   0.56  -0.17  -0.14  -0.45  -0.63   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.348, 0.5965, 0.7394, 0.8752, 0.9496, 0.9853, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.1    0.63  0.0   -0.46   0.1   -0.47  -0.39
 drpt          0.45   0.24  0.0    0.33  -0.77  -0.11  -0.19
 rcvr         -0.28  -0.54  0.0   -0.35  -0.44  -0.54   0.13
 trng          0.42   0.06  0.0   -0.73  -0.18   0.42   0.3
 drpt & trgt   0.58   0.0   0.0    0.12   0.35  -0.55   0.47
 rcvr & trgt  -0.45   0.49  0.0    0.13  -0.22   0.01   0.7
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + MixedModels.ZeroCorr((1 + a1 | Item))
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
4: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2007.1480                         
[2]         41  1980.2755  26.8726      13  0.0130
[3]         42  1969.2274  11.0480       1  0.0009
[4]         56  1954.9674  14.2600      14  0.4305
──────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0681020 0.2609636
drpt 0.0195454 0.1398050 -0.01
rcvr 0.0181474 0.1347124 -0.34 -0.45
trng 0.0299627 0.1730973 +0.21 +0.31 -0.11
drpt & trgt 0.0043532 0.0659788 -0.13 +0.53 -0.44 +0.43
rcvr & trgt 0.0019832 0.0445329 +0.56 -0.17 -0.14 -0.45 -0.63
trgt 0.0003244 0.0180111 . . . . . .
Item (Intercept) 0.0007211 0.0268527
a1 0.0004408 0.0209959 .
Residual 0.0724335 0.2691348

Perhaps not. We stay with LMM m_prm2

4.7 Compare models with goodness-of-fit statistics.

Code
let mods = [m_zcp2, m_prm2, m_prm1, m_cpx];
 DataFrame(;
    model=[:m_zcp2, :m_prm2, :m_prm1, :m_cpx],
    pars=dof.(mods),
    geomdof=round.(Int, (sum  leverage).(mods)),
    AIC=round.(Int, aic.(mods)),
    AICc=round.(Int, aicc.(mods)),
    BIC=round.(Int, bic.(mods)),
  )
end

4 rows × 6 columns

modelparsgeomdofAICAICcBIC
SymbolInt64Int64Int64Int64Int64
1m_zcp228293206320632252
2m_prm242279205320542336
3m_prm148274205920602382
4m_cpx56281206720682444

LMM m_prm2 is a defensible solution according to \(\Delta\) AIC, \(\Delta\) BIC suggests we should not bother with CPs.

Code
coeftable(m_prm2)
Coef. Std. Error z Pr(>
(Intercept) 7.69869 0.0371901 207.01 <1e-99
trgt -0.0118293 0.00451408 -2.62 0.0088
trng -0.070356 0.0262888 -2.68 0.0074
drpt 0.0457635 0.0231034 1.98 0.0476
rcvr -0.042934 0.0224154 -1.92 0.0554
a1 -0.108394 0.0326348 -3.32 0.0009
trgt & trng 0.00212125 0.00949196 0.22 0.8232
trgt & drpt -0.0370438 0.0150283 -2.46 0.0137
trgt & rcvr 0.00972961 0.0132374 0.74 0.4623
trgt & a1 -0.00701276 0.00399025 -1.76 0.0788
trng & a1 0.0279914 0.0230904 1.21 0.2254
drpt & a1 0.0273081 0.0201717 1.35 0.1758
rcvr & a1 -0.0167551 0.0195894 -0.86 0.3924
trgt & trng & a1 -0.00608894 0.00831708 -0.73 0.4641
trgt & drpt & a1 0.029943 0.0131747 2.27 0.0230
trgt & rcvr & a1 -0.00509192 0.011614 -0.44 0.6611

4.8 Diagnostic plots

Various visualizations are used to check whether or not data are defensibly modeled with an LMM. They may lead to removal of outliers, transformations of the dependent variable, and deliver valuable heuristic information to be followed up with exploratory post-hoc analyses or ideally replication of new insights gained this way. In practice, it appears that only severe violations will stop people from reporting a model.

4.8.1 Residuals over fitted

Code
scatter(fitted(m_prm2), residuals(m_prm2))

Looks like we missed some fast response times.

4.8.2 Q-Q plot

Code
qqnorm(m_prm2; qqline=:none)

Hm?

4.8.3 Residual distributions: observed vs. theoretical

Curves for residulas based on observed and theoretical values should correspond.

Code
let
  n = nrow(dat)
  dat_rz = (;
    value=vcat(residuals(m_prm2) ./ std(residuals(m_prm2)), randn(n)),
    curve=repeat(["residual", "normal"]; inner=n),
  )
  draw(
    data(dat_rz) *
    mapping(:value; color=:curve) *
    density(; bandwidth=0.1);
  )
end

Figure 1: Kernel density plot of the standardized residuals for model m1 versus a standard normal

They are a bit too narrow.

4.9 Conditional means of random effects

4.9.2 Borrowing-strength plots

Shrinkage refers to the adjustment of subject-level or item-level predictions by taking population estimates into account. The further a subject’s/item’s estimate is from the fixed effect or the more variable or less reliable the subject’s/item’s estimate, the more the prediction will be shrunk towards the population estimate. Alternative terms for shrinkage are “borrowing strength” (Tukey) and regularization. My favorite is actually Tukey’s because indeed we borrow strength from the population estimates to make predictions for individual subjects’ effects. The goal of this section to illustrate the results of borrowing strength.

Subject-related conditional means of random effects revealed information about individual differences beyond fixed effects. Would these results also be visible in unconditional means, that is when we compute GM and experimental effects within subjects (i.e., as fixed effects) without borrowing strength from the population estimates?

In the following plots, effect estimates based on alone on each subject’s data (i.e., no pooling of data, no borrowing of strength) are plotted in pink and the subjects’ conditional means shown in the caterpillar plots are plotted in blue. The arrows indicate how much a subject’s prediction is changed by borrowing strength from knowledge of the population estimates.

Code
shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm2)

Figure 3: Shrinkage plots of the subject random effects in model m_prm2

Outlier subject is gone.

Code
cm = raneftables(m_prm2);
sort(DataFrame(cm.Subj), ["(Intercept)", :trng, "drpt & trgt"])
sort(DataFrame(cm.Subj), :trng)

51 rows × 8 columns (omitted printing of 1 columns)

Subj(Intercept)drptrcvrtrngdrpt & trgtrcvr & trgt
StringFloat64Float64Float64Float64Float64Float64
1p14-0.5283230.159662-0.118234-0.3855480.0520919-0.00337878
2p760.064757-0.0954270.0965818-0.317894-0.0618810.0540763
3p360.373846-0.04461880.114011-0.277302-0.129930.109552
4p02-0.628903-0.1859060.239485-0.24184-0.0114483-0.0500986
5p390.292941-0.1233370.131867-0.189144-0.06628580.0443946
6p74-0.334411-0.0354860.058828-0.177365-0.00145978-0.0244133
7p670.01026580.00522903-0.0258083-0.1718380.005969190.0146027
8p590.3766580.0232177-0.253912-0.1379890.006891850.0472082
9p06-0.2388110.1408420.0320315-0.128487-0.04230110.0147919
10p540.319235-0.1623240.0453531-0.110458-0.0451120.0524742
11p42-0.1051710.295914-0.174955-0.1074080.0974469-0.00806991
12p23-0.0173974-0.07708850.064224-0.105342-0.09227020.0454142
13p30-0.34938-0.153186-0.0891651-0.0983450.00105422-0.0225567
14p45-0.374579-0.2942470.154624-0.089149-0.00555676-0.045438
15p84-0.432463-0.2158940.466653-0.0825766-0.170963-0.013087
16p83-0.23897-0.0537439-0.162418-0.07724240.0230343-0.0184981
17p720.177019-0.0778771-0.0337592-0.0503027-0.06298240.0679022
18p55-0.24107-0.0677181-0.0693261-0.0486936-0.0284252-0.0116131
19p27-0.013248-0.01112150.0152268-0.0443757-0.04412750.0253996
20p12-0.160657-0.04954520.0628843-0.0434829-0.0354820.00879487
21p690.0806806-0.0491015-0.0696128-0.03172180.0442580.025753
22p040.141053-0.04686010.0393693-0.0239661-0.01247610.0138384
23p770.150813-0.1014540.0587802-0.0141855-0.06364430.0327982
24p79-0.0709235-0.01864460.0306235-0.0134209-0.0220470.000434288
25p53-0.1129580.02493210.0212246-0.001522890.0380108-0.0234996
26p08-0.1603250.112630.03363740.005611310.0106376-0.000935882
27p250.2784660.141361-0.04764780.00933614-0.02496720.0311902
28p640.2660690.002937770.0343950.0113388-0.01736970.0242711
29p280.151958-0.0047019-0.1595940.01366860.0485822-0.00713976
30p570.00152955-0.0017727-0.01991410.01968260.037071-0.01744

5 Appendix

Code
versioninfo()
Julia Version 1.8.1
Commit afb6c60d69a (2022-09-06 15:09 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 8 on 8 virtual cores